/*
 * sm02.h
 *
 *  Created on: 2010.05.03
 *      Author: Tadas
 */

#ifndef SM02_H_
#define SM02_H_

void intPol(double *xa, double *ya, int n, double x, double &y, double &dy);

double simpson(double (*func)(double), double a, double b, int n){
    double x, h, h2, sum;

    int j;

    h = (b-a)/n;

    h2 = h/2;

    sum = func(a);

    for( j = 0; j < n; j++){
        x = a + j*h;
        sum += 4 * func(x+h2) + 2*func(x+h);
    }

    return (sum - func(b))*h2/3;
}

//s0 = [1/2*f(a) + 1/2*f(h)]h;
//s0 = [1/2*f(a) + 1/2*f(h/n)]*h/n;
//double trapecijos_tad( double (*func)(double), double a, double b, int n){
//    double sum, h;
//    h = (b - a)/n;
//    sum = 0;
//    for(int i = 0; i < n; i++){
//        //sum +=
//    }
//}

double trapec(double (*func)(double), double a, double b, int n){

    double x, tnm, sum, del;

    static double s;
    int j, i, it;

    //int it;

    if(n == 1){
        return (s = 0.5 * (b-a)*(func(a) - func(b)));
    }

    for(it = 1, i = 1; i < n - 1 and i < sizeof(it)*8; i++)
        it <<= 1;
    tnm = it;
    del = (b - a) / tnm;

    x = (a + 0.5*del);

    for(sum = 0.0, j = 0; j++ <= it; x += del){
        sum += func(x);
    }

    return (s = 0.5*s + sum*del);
}


double sumTrapec(double (*func)(double), double a, double b){
    int j, SNew, SOld = -1.0E-31;

    for( j = 1; j < JMAX; j++){
        SNew = trapec(func, a, b, j);

        if(fabs(SOld - SNew) < EPS*fabs(SOld)){
            return SNew;
        }

        if(SNew == 0.0 and SOld == 0.0 and j > 6){
            return 0.0;
            // return SNew;
        }

        SOld = SNew;
    }

    cout << "Per daug iteraciju; rezultato nera" << endl;
    return 0.0;
}

double sumtrapecSims(double (*func)(double), double a, double b){
    int j;

    double SOld, SNew, simNew, simOld;

    simOld = SOld = 1.0e-30;

    for(j = 1; j < JMAX; j++){
        simNew = trapec(func, a, b, j);

        SNew = (4.0 * simNew - simOld)/3.0;

        if(fabs(SNew - SOld) < EPS * fabs(SOld)){
            return SNew;
        }

        // tikrinam nulius

        SOld = SNew;

        simOld = simNew;
    }

    cout << "Nieko" << endl;

    return 0.0;
}


double sumTrapecRondberg(double (*func)(double), double a, double b){
    int j;
    double S, SErr;

    double hMas[JMAX], SMas[JMAX];

    hMas[0] = 1;

    for( j = 1; j < JMAX; j++){
        SMas[j - 1] = trapec(func, a, b, j);
        if(j >= K){
            intPol(hMas + j - K, SMas - j + K, K, 0.0, S, SErr);
            if(fabs(SErr) < EPS * fabs(S))
                return S;
        }

        hMas[j] = 0.25 * hMas[j - 1];
    }

    cout << "Daug zingsniu; nieko neturim" << endl;

    return 0.0;
}


/// Padarykite funkciją integruojančią
// Integruoti funkcijas knygutėje
// y = x^4 * log(x + sqrt(x^2 + 1))
// y = 10 * x, kai 0 < x < 0.5
//   = (20 - 10*x)/3, kai 0.5 <= x < 2
// palyginti našumą

#endif /* SM02_H_ */
